home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-04
/
extenpr.zip
/
EXTEND.FOR
< prev
next >
Wrap
Text File
|
1988-06-04
|
16KB
|
709 lines
C
C This is a collection of subroutines that perform extended precision
C floating point arithmetic.
C
C By Mark P. Esplin
C 16 Shawsheen Rd.
C Billerica, Ma 01821
C
C Whole programs do not have to be written using extended
C precision. Only the part that needs extra precision need use it,
C then convert back to regular REAL*8 numbers.
C The routines for addition and multiplication use the same algorithms
C that are taught in grade school but using base 10000 numbers.
C This is not the most efficient way to do it, but
C it makes it possible to use standard FORTRAN77 (almost).
C These routines should work with any computer and compiler that supports
C INTEGER*2 and INTEGER*4 arithmetic.
C Each number is stored in an INTEGER*2 array of length N.
C The first word, IA(1), is the exponent in base 10000 digits.
C The least significant digit is in IA(2) and the most significant
C digit is in IA(N). If IA(N) > 5000 the number is negative.
C When IA(N) = 5000 the number is assumed invalid.
C Base 10000 digits have worse roundoff characteristics than binary.
C Since roundoff causes the last base 10000 digit to be inaccurate,
C that really means the last 4 decimal digits are inaccurate. Be
C sure to use more digits than is necessary to allow for roundoff.
C N must be at least 6 or some routines won't work correctly.
C
C Subroutines with their FORTRAN equivalent in this package are:
C
C ADD(A,B,C,N) C=A+B
C NEG(A,N) A=-A
C MUL(A,B,C,N) C=A*B
C INV(A,B,N,W) B=1/A W is work area of size 3*N
C EXSQRT(A,B,N,W) B=SQRT(A) W is work area of size 4*N
C EXHALF(A,N) A=1/2 much faster than using INV
C PRTEXT(A,N) WRITE(*,*) A
C WRTEXT(IU,A,N) WRITE(IU,*) A
C CVI2EX(I2,A,N) A=I2 I2 is INTEGER*2
C CVR8EX(R,A,N) A=R R is REAL*8
C CVEXR8(A,N,R) R=A R is REAL*8
C COPY(A,B,N) B=A
C IEXCMP(A,B,N,W) Compares two positive numbers A and B returns 1, 0, -1
C CHGPRE(A,NA,B,NB) Changes the precision of A form NA to NB
C
C There are other internal routines, but they would probably be of limited
C value for most users.
C
C
SUBROUTINE INCEX(IA,N)
INTEGER*2 IA(*),N
C This routine increments the mantissa of extended number IA.
C This routine is used by other routines to perform rounding.
INTEGER*2 IT,ICAR,I
ICAR=1
DO 10 I=2,N
IT=IA(I)+ICAR
ICAR=IT/10000
IA(I)=IT-ICAR*10000
C In almost all cases this routine will end after a few passes.
IF(ICAR.EQ.0) GOTO 20
10 CONTINUE
C Shift back if necessary
C The only time that a shift is necessary is when IA(N) = 5000 and
C the rest of the mantissa is zero.
20 CONTINUE
IF(IA(N).EQ.5000) THEN
IA(1)=IA(1)+1
DO 30 I=2,N-1
IA(I)=IA(I+1)
30 CONTINUE
IA(N)=0
C Don't need to round again since if a shift was necessary
C the least significant digit is zero.
ENDIF
RETURN
END
C
C
C
C
C
FUNCTION ICKCON(IA,N)
C This routine is used to check if an extended number IA is converging
C to an integer value.
C This routine would probably not be very useful for the average user.
C The value returned is 1 if IA has converged. The list significant digit
C is not checked since it might be just round off.
C This function return a 1 in IA has converged to an integer otherwise
C it returns a 0.
INTEGER*2 IA(*)
INTEGER*2 N,I
INTEGER*2 ITEST
ICKCON=0
C The number must have already converged past the N-2 th digit before
C this routine is called.
ITEST=IA(N-2)
DO 10 I=3,N-3
IF(IA(I).NE.ITEST) RETURN
10 CONTINUE
ICKCON=1
RETURN
END
SUBROUTINE NEG(IA,N)
C This routine changes the sign of extended precision number IA.
INTEGER*2 IA(*),N
INTEGER*2 I,J
J=2
DO 10 I=2,N
IF(IA(J).NE.0) GOTO 20
J=J+1
10 CONTINUE
C Number was zero
RETURN
C The first non-zero digit has been reached
20 IA(J)=10000-IA(J)
J=J+1
DO 30 I=J,N
IA(I)=9999-IA(I)
30 CONTINUE
RETURN
END
C
C
C
C
SUBROUTINE PRTEXT(IA,N)
C This routine writes extended precision numbers to output.
C This routine could use some work to make a more readable output.
INTEGER*2 IA(*),N
INTEGER*2 I,INEG
INEG=0
WRITE(*,*) ' '
IF(IA(N).GT.5000) THEN
WRITE(*,*) 'Minus'
CALL NEG(IA,N)
INEG=1
ENDIF
WRITE(*,900) IA(1)
900 FORMAT(' Exponent ',I6)
WRITE(*,910) (IA(I),I=N,2,-1)
910 FORMAT(10I5)
C Restore IA to previous value.
IF(INEG.NE.0) CALL NEG(IA,N)
RETURN
END
C
C
C
C
SUBROUTINE WRTEXT(IUNIT,IA,N)
C This routine writes extended precision numbers to file IUNIT.
C This routine could use some work to make a more readable output.
INTEGER*2 IUNIT,IA(*),N
INTEGER*2 I,INEG
INEG=0
WRITE(IUNIT,*) ' '
IF(IA(N).GT.5000) THEN
WRITE(IUNIT,*) 'Minus'
CALL NEG(IA,N)
INEG=1
ENDIF
WRITE(IUNIT,900) IA(1)
900 FORMAT(' Exponent ',I6)
WRITE(IUNIT,910) (IA(I),I=N,2,-1)
910 FORMAT(10I5)
C Restore IA to previous value.
IF(INEG.NE.0) CALL NEG(IA,N)
RETURN
END
C
C
C
C
C
SUBROUTINE NORMIZ(IA,N)
C This routine is used by other routines to normalize
C floating point number after addition and other such processes.
INTEGER*2 IA(*),N
INTEGER*2 I,ISHFT
IF(IA(N).GT.5000) THEN
C This is a negative number.
DO 10 I=N,2,-1
IF(IA(I).NE.9999) GOTO 20
10 CONTINUE
I=I+1
20 ISHFT=N-I
IF(IA(I).LE.5000) ISHFT=ISHFT-1
ELSE
C This is a positive number.
DO 30 I=N,2,-1
IF(IA(I).NE.0) GOTO 40
30 CONTINUE
C This number is zero
RETURN
40 ISHFT=N-I
IF(IA(I).GE.5000) ISHFT=ISHFT-1
ENDIF
IF(ISHFT.EQ.0) RETURN
IA(1)=IA(1)-ISHFT
C Shift left the appropriate amount.
DO 50 I=N,ISHFT+2,-1
IA(I)=IA(I-ISHFT)
50 CONTINUE
C Zero unknown digits.
DO 60 I=2,ISHFT+1
IA(I)=0
60 CONTINUE
RETURN
END
C
C
C
C
C
SUBROUTINE ADD(IA,IB,IC,N)
C This routine performs extended precision addition.
C IC = IA + IB. IA and IB are not changed.
INTEGER*2 IA(*),IB(*),IC(*)
INTEGER*2 N,IPA,IPB,NML,I
INTEGER*2 ISGNA,ISGNB,ISGNC,ISGEXT
INTEGER*2 IT,ICAR
ISGNA=1
IF(IA(N).GT.5000) ISGNA=-1
ISGNB=1
IF(IB(N).GT.5000) ISGNB=-1
C Figure out pointers for main loop
IF(IA(1).GE.IB(1)) THEN
IC(1)=IA(1)
IPA=2
IPB=IA(1)-IB(1)+2
NML=N-IPB+2
IF(NML.LT.2) THEN
C The difference between the two numbers is too large.
CALL COPY(IA,IC,N)
RETURN
ENDIF
IF(IPB.EQ.2) THEN
ICAR=0
ELSE
ICAR=(IB(IPB-1)+5000)/10000
ENDIF
ELSE
IC(1)=IB(1)
IPB=2
IPA=IB(1)-IA(1)+2
NML=N-IPA+2
IF(NML.LT.2) THEN
C The difference between the two numbers is too large.
CALL COPY(IB,IC,N)
RETURN
ENDIF
ICAR=(IA(IPA-1)+5000)/10000
ENDIF
C Main addition loop.
DO 10 I=2,NML
IT=IA(IPA)
IT=IT+IB(IPB)+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
IPA=IPA+1
IPB=IPB+1
10 CONTINUE
C Do sign extended part.
IF(IA(1).GT.IB(1)) THEN
IF(ISGNB.GE.1) THEN
ISGEXT=0
ELSE
ISGEXT=9999
ENDIF
DO 20 I=NML+1,N
IT=IA(I)
IT=IT+ISGEXT+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
20 CONTINUE
ELSE
IF(ISGNA.GE.1) THEN
ISGEXT=0
ELSE
ISGEXT=9999
ENDIF
DO 30 I=NML+1,N
IT=IB(I)
IT=IT+ISGEXT+ICAR
ICAR=IT/10000
IC(I)=IT-ICAR*10000
30 CONTINUE
ENDIF
C Shift right one place if it has overflowed during addition.
IF(ISGNA*ISGNB.EQ.1) THEN
ISGNC=1
IF(IC(N).GE.5000) ISGNC=-1
IF(ISGNA*ISGNC.EQ.-1) THEN
C An overflow has occurred.
IC(1)=IC(1)+1
ICAR=IC(2)
DO 40 I=2,N-1
IC(I)=IC(I+1)
40 CONTINUE
IC(N)=ISGEXT
C Do rounding and return.
IF(ICAR.GE.5000) CALL INCEX(IC,N)
RETURN
ENDIF
ENDIF
C Normalize answer.
CALL NORMIZ(IC,N)
RETURN
END
C
C
C
C
SUBROUTINE MUL(IA,IB,IC,N)
C This routine multiplies two extended precision numbers.
C N must be at least 4 or this routine won't work.
C IC = IA*IB
INTEGER*2 IA(*),IB(*),IC(*),N
INTEGER*2 ISNA,ISNB,IXC1,IXC2,ISHFT
C The computation is carried out with two extra placesIXC1 and IXC2
C This makes it so precision is not lost if the answer has to be
C shifted left. IXC1 is the least significant.
INTEGER*4 IT,ICAR
INTEGER*2 I,J,K,IAS
ISNA=1
ISNB=1
IF(IA(N).GT.5000) THEN
CALL NEG(IA,N)
ISNA=-1
ENDIF
IF(IB(N).GT.5000) THEN
CALL NEG(IB,N)
ISNB=-1
ENDIF
IC(1)=IA(1)+IB(1)
IXC1=0
IXC2=0
IAS=N+1
C Step through each IB
DO 30 I=2,N
ICAR=0
IF(I.LE.(N-2)) THEN
C Round digit that is past least significant.
IT=IA(IAS-3)
IT=IT*IB(I)+5000
ICAR=IT/10000
ENDIF
IF(I.LE.(N-1)) THEN
C Do extended precision 1
IT=IA(IAS-2)
IT=IT*IB(I)+ICAR+IXC1
ICAR=IT/10000
IXC1=IT-ICAR*10000
ENDIF
C Do extended precision 2
IT=IA(IAS-1)
IT=IT*IB(I)+ICAR+IXC2
ICAR=IT/10000
IXC2=IT-ICAR*10000
C The first time through the I loop the J loop is not executed
K=2
DO 20 J=IAS,N
IT=IA(J)
IT=IT*IB(I)+ICAR+IC(K)
ICAR=IT/10000
IC(K)=IT-ICAR*10000
K=K+1
20 CONTINUE
IC(K)=ICAR
IAS=IAS-1
30 CONTINUE
C Normalize output if needed rounding when possible.
IF((IC(N).EQ.0).AND.(IC(N-1).LT.5000)) THEN
IF((IC(N-1).EQ.0).AND.(IC(N-2).LT.5000)) THEN
ISHFT=2
ELSE
ISHFT=1
ENDIF
IC(1)=IC(1)-ISHFT
DO 40 I=N,ISHFT+2,-1
IC(I)=IC(I-ISHFT)
40 CONTINUE
IF(ISHFT.EQ.2) THEN
C No rounding is possible in this case.
IC(3)=IXC2
IC(2)=IXC1
ELSE
IC(2)=IXC2
IF(IXC1.GE.5000) CALL INCEX(IC,N)
ENDIF
ELSE
IF(IXC2.GE.5000) CALL INCEX(IC,N)
ENDIF
C Set sign
IF(ISNA*ISNB.LT.0) CALL NEG(IC,N)
RETURN
END
C
C
C
C
SUBROUTINE EXSQRT(IA,IB,N,IW)
C This routine does an extended precision square root.
C IB = SQRT(IA)
C It uses an iterative method. See Arithmetic Operations in Digital
C Computers by R. K. Richards D. Van Nostrand Company, Inc. Copyright 1955.
C B(k+1) = B(k)/2 * (3 - B(k)^2/A) This series only requires one division.
C IW is a work area of size 4*N
INTEGER*2 IA(*),IB(*),IW
INTEGER*2 N,NUSED,NUSEDL
INTEGER*2 IALD,IMONE,ITHREE
DIMENSION IW(N,4)
REAL*8 R
IALD=0
IMONE=-1
ITHREE=3
IF(IA(N).GE.5000) THEN
WRITE(*,*) 'Square root of negative number?'
STOP
ENDIF
C Find first guess.
CALL CVEXR8(IA,N,R)
R=SQRT(R)
NUSED=6
CALL CVR8EX(R,IB,NUSED)
C IW(1) is 1/IA at full precision.
CALL INV(IA,IW(1,1),N,IW(1,2))
C Beginning of iterative loop
10 NUSEDL=NUSED
NUSED=(NUSED-1)*2
IF(NUSED.GT.N) THEN
NUSED=N
IALD=1
ENDIF
CALL CHGPRE(IB,NUSEDL,IB,NUSED)
CALL MUL(IB,IB,IW(1,2),NUSED)
C Return if convergence is complete
IF(IALD.EQ.1) THEN
IF(IEXCMP(IA,IW(1,2),N,IW(1,3)).EQ.0) RETURN
ENDIF
CALL CHGPRE(IW(1,1),N,IW(1,3),NUSED)
CALL MUL(IW(1,2),IW(1,3),IW(1,4),NUSED)
CALL NEG(IW(1,4),NUSED)
CALL CVI2EX(ITHREE,IW(1,2),NUSED)
CALL ADD(IW(1,2),IW(1,4),IW(1,3),NUSED)
C IW(3) is now equal to 3 - B(k)^2/A
CALL EXHALF(IW(1,2),NUSED)
CALL MUL(IW(1,2),IB,IW(1,4),NUSED)
CALL MUL(IW(1,4),IW(1,3),IB,NUSED)
GOTO 10
END
C
C
C
C
C
C
SUBROUTINE CVI2EX(I2,IA,N)
C This routine converts an INTEGER*2 number to extended precision.
INTEGER*2 I2,IA(*)
INTEGER*2 N,I
INTEGER*2 ISGN,ITEMP
C Zero array before start
DO 10 I=1,N
IA(I)=0
10 CONTINUE
ITEMP=I2
ISGN=1
IF(ITEMP.LT.0) THEN
ISGN=-1
ITEMP=-ITEMP
ENDIF
IF(ITEMP.GE.5000) THEN
IA(1)=2
IA(N)=ITEMP/10000
IA(N-1)=ITEMP-IA(N)*10000
ELSE
IA(1)=1
IA(N)=ITEMP
ENDIF
IF(ISGN.EQ.-1) CALL NEG(IA,N)
RETURN
END
C
C
SUBROUTINE COPY(IA,IB,N)
C This routine copies extended number IA to IB.
INTEGER*2 IA(*),IB(*),N
INTEGER*2 I
DO 10 I=1,N
IB(I)=IA(I)
10 CONTINUE
END
C
C
C
C
SUBROUTINE EXHALF(IA,N)
C This routine returns an extended precision 1/2.
INTEGER*2 IA(*)
INTEGER*2 N,I
IA(1)=1
IA(N)=0
IA(N-1)=5000
DO 10 I=2,N-2
IA(I)=0
10 CONTINUE
RETURN
END
C
C
C
C
C
SUBROUTINE CVR8EX(R,IA,N)
C Convert REAL*8 R to an extended precision number IA.
C N must be greater than or equal to 6.
REAL*8 R,RT
INTEGER*2 IA(*)
INTEGER*2 N,I
INTEGER*2 ISIGN
ISIGN=1
RT=R
IF(R.LT.0) THEN
ISIGN=-1
RT=-RT
ENDIF
IA(1)=LOG(RT)/LOG(10000.)
RT=RT/(10000.**IA(1))
IA(1)=IA(1)+1
IF(RT.GE.5000) THEN
IA(1)=IA(1)+1
RT=RT/10000.
ENDIF
DO 10 I=N,N-4,-1
IA(I)=RT
RT=RT-IA(I)
RT=RT*10000
10 CONTINUE
DO 20 I=N-5,2,-1
IA(I)=0
20 CONTINUE
IF(ISIGN.LT.0) CALL NEG(IA,N)
RETURN
END
C
C
C
SUBROUTINE CVEXR8(IA,N,R)
C Convert an extended precision number IA to REAL*8 R.
C N must be greater than or equal to 6.
REAL*8 R,RT,EXP
INTEGER*2 IA(*)
INTEGER*2 N,I
INTEGER*2 ISIGN
ISIGN=1
IF(IA(N).GE.5000) THEN
ISIGN=-1
CALL NEG(IA,N)
ENDIF
RT=0
EXP=1
EXP=(EXP/10000)**5
DO 10 I=N-4,N
RT=RT+IA(I)*EXP
EXP=EXP*10000
10 CONTINUE
EXP=10000
R=RT*(EXP**IA(1))
IF(ISIGN.LE.0) R=-R
RETURN
END
C
C
C
C
SUBROUTINE CHGPRE(IA,NA,IB,NB)
C This routine changes the precision of an extended precision number.
C This routine is useful for many iterative solutions where the full
C precision is not needed until the final iteration.
C IA and IB can be the same array.
C Input is IA of size NA, output is IB of size NB.
INTEGER*2 IA(*),IB(*),ICAR
INTEGER*2 NA,NB,I,J
IB(1)=IA(1)
IF(NA.GT.NB) THEN
C IA has the most precision.
ICAR=0
J=NA-NB+1
IF(IA(J).GT.5000) ICAR=1
DO 10 I=2,NB
J=J+1
IB(I)=IA(J)
10 CONTINUE
IF(ICAR.EQ.1) CALL INCEX(IB,NB)
ELSE
C IB has the most precision.
J=NB
DO 20 I=NA,2,-1
IB(J)=IA(I)
J=J-1
20 CONTINUE
DO 30 I=2,J
IB(I)=0
30 CONTINUE
ENDIF
RETURN
END
C
C
C
C
C
SUBROUTINE INV(IX,IY,N,IW)
C This routine does an extended precision inverse.
C IY=1/IX
C It uses an iterative method. A routine more like MUL would probably
C be more efficient and would certainly take less storage, but I haven't
C gotten around to writing it yet.
C Y(I+1)=Y(I)*(2-X*Y(I))
C IW is a work area of size 3*N
INTEGER*2 IX(*),IY(*),IW
INTEGER*2 N,NUSED,NUSEDL
INTEGER*2 IALD,ISGN,INEGT
DIMENSION IW(N,3)
REAL*8 R
INEGT=-2
IALD=0
ISGN=1
IF(IX(N).GE.5000) THEN
CALL NEG(IX,N)
ISGN=-1
ENDIF
C Find first guess.
CALL CVEXR8(IX,N,R)
R=1/R
NUSED=6
CALL CVR8EX(R,IY,NUSED)
10 NUSEDL=NUSED
NUSED=(NUSED-1)*2
IF(NUSED.GT.N) THEN
NUSED=N
IALD=1
ENDIF
CALL CHGPRE(IX,N,IW(1,1),NUSED)
CALL CHGPRE(IY,NUSEDL,IY,NUSED)
CALL MUL(IW(1,1),IY,IW(1,2),NUSED)
C Return if convergence is complete
IF(IALD.EQ.1) THEN
IF(ICKCON(IW(1,2),N).EQ.1) THEN
IF(ISGN.LT.0) CALL NEG(IY,N)
RETURN
ENDIF
ENDIF
CALL CVI2EX(INEGT,IW(1,1),NUSED)
CALL ADD(IW(1,1),IW(1,2),IW(1,3),NUSED)
CALL NEG(IW(1,3),NUSED)
CALL MUL(IY,IW(1,3),IW(1,1),NUSED)
CALL COPY(IW(1,1),IY,NUSED)
GOTO 10
END
C
C
C
C
C
C
C
FUNCTION IEXCMP(IA,IB,N,IW)
C This routine compares two extended precision numbers.
C The returned value is 1 if IA > IB, 0 if IA is within roundoff of IB, etc.
C IW is a work area of size N.
C Only works for positive numbers
INTEGER*2 IA(*),IB(*),IW(*)
INTEGER*2 N,I
REAL*8 A,B
REAL DIF
C First check if the numbers are close.
CALL CVEXR8(IA,N,A)
CALL CVEXR8(IB,N,B)
DIF=(A-B)/(A+B)
IF(ABS(DIF).GT.1E-12) THEN
IEXCMP=1
IF(DIF.LT.0) IEXCMP=-1
RETURN
ENDIF
C The values are close
CALL NEG(IB,N)
CALL ADD(IA,IB,IW,N)
CALL NEG(IB,N)
C Check for zero
IF(IA(1)-IW(1).GE.N-2) THEN
C Within roundoff of zero
IEXCMP=0
RETURN
ENDIF
DO 10 I=2,N
IF(IW(I).NE.0) GOTO 20
10 CONTINUE
IEXCMP=0
RETURN
C Results are not zero
20 IEXCMP=1
IF(IW(N).GT.5000) IEXCMP=-1
RETURN
END